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Abstract 

A great improvement to the insight on brain function that we can get from 
fMRI data can come from effective connectivity analysis, in which the flow of 
information between even remote brain regions is inferred by the parameters 
of a predictive dynamical model. As opposed to biologically inspired mod- 
els, some techniques as Granger causality (GC) are purely data-driven and 
rely on statistical prediction and temporal precedence. While powerful and 
widely applicable, this approach could suffer from two main limitations when 
applied to BOLD fMRI data: confounding effect of hemodynamic response 
function (HRF) and conditioning to a large number of variables in presence 
of short time series. For task-related fMRI, neural population dynamics can 
be captured by modeling signal dynamics with explicit exogenous inputs; for 
resting-state fMRI on the other hand, the absence of explicit inputs makes 
this task more difficult, unless relying on some specific prior physiological 
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hypothesis. In order to overcome these issues and to allow a more general 
approach, here we present a simple and novel blind-deconvolution technique 
for BOLD-fMRI signal. In a recent study it has been proposed that relevant 
information in resting-state fMRI can be obtained by inspecting the discrete 
events resulting in relatively large amplitude BOLD signal peaks. Follow- 
ing this idea, we consider resting fMRI as 'spontaneous event-related', we 
individuate point processes corresponding to signal fluctuations with a given 
signature, extract a region-specific HRF and use it in deconvolution, after 
following an alignment procedure. Coming to the second limitation, a fully 
multivariate conditioning with short and noisy data leads to computational 
problems due to overfitting. Furthermore, conceptual issues arise in pres- 
ence of redundancy. We thus apply partial conditioning to a limited subset 
of variables in the framework of information theory, as recently proposed. 
Mixing these two improvements we compare the differences between BOLD 
and deconvolved BOLD level effective networks and draw some conclusions. 

Keywords: BOLD signal, Deconvolution, Effective connectivity, Granger 
causality 



1. Introduction 



We can learn a lot on the functioning of the human brain in health and 
disease when we consider it as a large-scale complex n etwork, whose proper- 
ties c an be analyzed using graph theoretical analysis (IBullmore and Spornsl . 
20091 ). With the advent of miscellaneous and noninvasive MRI techniques, 
this connectome has been mainly characterized by either structural or func- 
tional connectivity. Structural connectivity is c ommonly based on whi te mat- 
ter tracts quantified by diffusion tractography ( Hagmann et all 2008 ); func- 
tional connectivity relie s on the other hand o n statistical dependencies such 
as temporal correlation ([Salvador et all 120051 ) . An impo rtant addition to this 
fram ework can come from effective connectivity analysis fjStephan and Roebroeck 
20121 ). in which the flow of information between even remote brain regions is 
inferred by the parameters of a predictive dynamical model. 



For some techniques, su ch as dynamic causal model l ing (PCM) an d struc - 



tural equation modelling (jBuchel and Fristonl . Il997t iFriston et all 120031 ) . 
these models are built and validated from specific anatomical and physiolog- 
ical hypotheses. Other techniques such as Granger causality analysis (GCA) 
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(IBressler and Sethi . 1201 ll ). are on the other hand data-driven and rely purely 
on statistical prediction and temporal precedence. While powerful and widely 
applicable, this last approach could suffer from two main limitations when ap- 
plied to blood-oxygenation level-dependent (BOLD)-functional MRI (fMRI) 
data: confounding effect of hemodynamic response function (HRF) and con- 
ditioning to a large number of variables in presence of short time series. 
Early interpretation of fMRI based directed connectivity by GCA always as- 
sumed homogeneous hemodynamic processes over the brain; several studies 
have pointed out that this is indeed not the case and that we are faced with 
variab le HRF latency across physiological processes a nd distinct brain re- 
gions ( iRoebroeck et all 1201 ll ; IValdes-Sosa et all 1201 if ). Recently, a number 
of studies have ad dressed this issue proposing to model the HRF according 



to several recipes (IBakhtiari and Hossein-Zadehl . |2012| ; lHavlicek et all 12011 



20101 ; iRyali et all 1201 if ). As well, a recent study has proposed that it would 



still feasible to infer connectivity at BOLD level, under t he assumption that 
Gran ger causality is theoretically invariant under filtering f lBarnett and Sethi . 
20 111 ) and that the HRF can be considered as a filter. It is still unclear 
whether and how specific effects related to HRF disturb the inference of 
temporal precedence. In addition a simulated or experimental ground truth 
is difficult to obtain, though some studies on simulated fMRI data have 
tried to re veal the relationship be t ween neural- l evel a nd BOLD-level causal 



influence (IDeshpande et all |2010| ; ISmith et all 1201 if ). A considerable help 



to obtain the HRF for deconvolution could come from multimodal imaging 
where the high temporal resolution of EEG is combined to the high spatial 
resolution of fMRI, but this experimental approach is still far from being 
optimal and w idely applicable. HRF ha s been studied almost since the early 



days of fMRI (IHandwerker et all |2012j). For task- related fMRI, neural pop 



ulation dynamics can be captured by modeling signal dyn amics with explicit 



exogenous inputs ( jFriston et all . l2008t iRiera et all 120041 1. i.e. deconvolution 



according to the explicit task design is possible in this case (IBuxton et al. 



1998; Friston et al. 



20001 ; iGloverl . Il999l ). For resting-state fMRI on the other 



hand, the absence of explicit inputs makes this task more difficult, unless 



relying on some spec ific prior physiological hypothesis (jFriston et all 12008 



Havlicek et all 1201 if ) . In order to overcome these issues and to allow a more 
general approach, here we present a simple and novel blind-deconvolution 
technique for BOLD-fMRI signal. 



Coming to the second limitation, in order to distinguish among direct 
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and mediated influences in multivariate datasets it is necessary to condi- 
tion the analysis to other variables. A bivariate analysis would indeed lead 
to the detection of many false positives. In presence of a large number of 
variable and short time series, a fully multivariate conditioning could lead 
to computational problems due to the overfitting. Fur thermore, conceptual 



issues would arise in p resence of redundant variables flAngelini et al.l . 12010 



Marinazzo et all 1201 0l ). In this paper we thus apply partial conditioning for 
Grang er Causality iPCGCfl to a limited subset of variables, as recently pro- 
posed ( iMarinazzo et all 120121 ) for reconstructing the BOLD and deconvolved 
BOLD level effective connectivity network (ECN) and compare them. 



2. Materials and methods 



2.1. Blind- deconvolution in resting-state fMRI data 



(David et al.. 


2008: 


Glover. 


1999) 



formation from neural activation to BOLD response can be modeled as a 
linear and time invariant system, measured fMRI data b(t) can be seen as 
the result of the convolution of neural states s(t) with a HRF h(t): 



b(t) = s(t) ® h(t) + e(i) 



(1) 



Where t is the time and ® denotes convolution. e(t) is the noise in the 
measurement, which we assume to be white. Since the right side of the 
above equation includes three unobservable quantities, in order to solve the 
equation for h{t) we need to substitute s(t) with a hypothetical model of the 
neural activation for s(t). Here we employ a simple on-off model of activation 
to model s(t): 



(2) 



where 8(t — r) is the delta function. This allows to fit the HRF h(t) according 
to s(t) using a canonical HRF (two gamma functions) and two derivatives 
(mult ivariate Taylo r expa nsion: temporal derivative and dispersion deriva- 
tive) (IFriston et all 120001 ) . as is common in most fMRI studies. 



1 Plcase note that this approach is different from partial Granger causality (PGC) (Guo 
et al. (2008), Journal of Neuroscience Methods, 172, 79) 
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Once calculated h(t), we can obtain an approximation s(t) of the neural 
signal from the observed data using a Wiener filter 



s(t) = d(t) ® b(t) 



(3) 



Let H(oj), B(u), E(u), and D(u) be the Fourier transforms of h(t), b(t), e(t), 
and d(t), respectively. Then 



D(u) 



H*(u) 



(4) 



where * denotes complex conjugate. The estimation s(t) of the neural states 
s(t) is then given by 



s(t) = FT' 1 {D(co)B(co)} = FT' 1 j 



H*(u)B(u}) 



H(lo)\ 2 + \E(uj)\ 2 



(5) 



Where FT 1 is the inverse Fourier transform operator. 



For task-related fMRI, the stimulus function provides the prior expec- 
tations about neural activity and a generative model whose inversion cor- 
responds to deconvolution; this is in principle not the case for resting-state 
fMRI. Nonetheless there is increasing evidence of spe cific events and neural 



states that govern th e dynamics of the brain at rest (IDeco and Jirsal . 12012 



Petridou et all . 120121 ) . Furthermore, Tagliazucchi et al. proposed that these 
events are reflected by relatively large amplitude BOLD signal peaks and 
thus that such fluctuat ions could encode releva nt information from resting- 
state fMRI recordings fjTagliazucchi et all 120121 ) . Inspired by their work, we 
consider resting-state fMRI as spontaneous event-related, and we propose to 
extract the HRF from those pseudo-events. After doing this, we can employ 
the deconvolution model in the same way as described above. It is known 
that the BOLD response is much slower than the neural activation that is 
presumed to drive it. Consequently, the peak of the BOLD signal lags behind 
the peak of neural activation (i.e. by k points). So here we assume that these 
events are generated from s(t). 

Glover pointed out that the noise spectrum in task-related fMRI can 
be obta i ned fr om time series measurements in nonactivated cortical regions 
( IGloverl . 119991 ); here we extend the model to cope with resting-state fMRI for 
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which there is no explicit activation. In this study we assumed covariance of 
noise e equal to cov [b(t) — s(t) <g) h(t)]. 

In order to obtain a value for k, we search all integer values in the interval 
[0 Umax], where K max is an arbitrary maximum value, choosing the one for 
which the noise error covariance is smallest as the onset. By this method we 
can perform deconvolution on all BOLD signals, requiring no information on 
timing or a priori spatial information of events; furthermore, the time series 
could be the average of time series over a region of interest with any scale, 
or series extracted by independent or principal component analysis. A flow 
chart for BOLD signal deconvolution is shown in FigJTJ 
This is the pseudo-code for our procedure. 

Pseudo-code 

i Preprocess time series (e.g. detrend, normalize etc.). 

ii Find a time set S in which the BOLD values exceed a given 
threshold around a local maximum. 

iii choose a maximum time delay K max 
FOR n = to K max 

S n = S — n 

s n (t) = l,te S n ; s n (t) = 0, t S n . 

Fit a general linear model using s n and canonical HRF with 
time and dispersion derivatives. 

END FOR 

iv Let e K = min {e n }, where e n = cov [b(t) — s n (t) <g> h n (t)]. 

0^71^. K, m ax 

v Follow equation 4 and 5, using HRF h K and e K for deconvolution, 
get s(t). 
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Figure 1: Flow chart for blind-deconvolution procedure. 1. the pre-processed (detrended 
and normalized) observed BOLD signal is evaluated against a given threshold obtaining 
several sets of putative onsets for pseudo-events. 2. the time deviation of the timing 
sets is adjusted; the set with smallest noise error covariance will represent the event. 3. 
the observed BOLD signal is deconvolved into a neural signal by using the corresponding 
HRF. 



2.2. Partially conditioned Granger causality 

Here we employ a methodology proposed in (IMarinazzo et all 120121 ) which 
allows to compute Granger causality conditioned to a limited number of vari- 
ables in the framework of information theory. The idea is that conditioning 
on a small number of variables, chosen as the most informative for the can- 
didate driver variable, is sufficient to remove indirect interactions for sparse 
connectivity patterns. 

We consider n covariance-stationary variables {xi(t)} i=1 ... n , denoting the 
state vectors as: 

X a {t) = (x a {t-m) r -- ,x a (t-l)) (6) 

m being the model order. Let e(x a \Y) be the mean squared error prediction of 
x a on the basis of the vectors Y. The partially conditioned Granger causality 
index c{fi — > a) is defined as follows: 

C(/3 a) = l ° 9 e(x a \ZUX,) (?) 
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Where Z = {X ix , • ■ ■ , X in } is a set of the variables, in {X±, X 2 , ■ ■ ■ , X n }\Xf- 
which are most informative for Xp. We adopt the following approximate 
strategy for Z : given the previous Z k _i, the set Z k is obtained adding the 
variable with greatest information gain. This is repeated until variables 
are selected. 

2.3. Simulation Datasets: NetSim 

A method for establishing a ground truth for fMRI data ha s not reached 



a gen eral consensus. Recently a benchmark dataset, NetSim (ISmith et al. 



20 111 ) has attracted a lot of attention. Previous studies have shown that lag- 
based methods perform very poorly on these datasets; it is anyway worthy to 
mention that these data are simulated under the DCM framework, contain no 
reciprocal connections and only Gaussian noise, limiting their universality as 
ground truth. Here we analyzed the largest of these datasets, consisting of 50 
nodes. After deconvolution the sensitivity improved significantly, increasing 
from 20% to 30%. Also the specificity improved from 88% to 94%. This 
does not render GC the method of choice for these data, for which we also 
have to point out that "neural events" and noise are not distinguishable, but 
gives nonetheless an indicative result for the usefulness of deconvolution of 
the BOLD signal. 

2.4- Resting-State fMRI Datasets 

In order to investigate the role of repetition time (TR) on the deconvo- 
lution procedure and on the effective network reconstruction, our analyses 
were performed on a resting-state fMRI dataset which has been publicly re- 
leased in the "1000 Functional Connectomes Project"!!. All participants had 
no history of neurological and psychiatric disorders and all gave the informed 
consent approved by local Institutional Review Board. During the scanning 
participants were instructed to keep their eyes closed, not to think of any- 
thing in particular, and to avoid falling asleep. 

Two data sets with different TR (TR=1.4s and TR=0.645s) were acquired 
on Siemens 3T Trio Tim scanners using developed multiplexed echo planar 



2 http://fcon_1000. projects. nitre. org, accessed march 2012. 
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imaging (IFeinberg et al.l . |2010| ). As specified in detail below, two resting- 
state fMRI data are included in the protocol - a TR=0.645s (3mm isotropic 
voxels, 10 minutes) to provide optimal temporal resolution and TR=1.4s 
(2mm isotropic voxels, 10 minutes) to provide optimal spatial resolution. 
The third data set, acquired on a 4T scanner, contains standard resting- 
state fMRI acquisitions with a longer TR (TR=3, 4mm isotropic voxels, 
5 minutes). For more detail on subject and data information, please see 
websiteP @. 



2. 5. Data preprocessing 

Preprocessing of resting-state images was performed using the Statisti- 



cal Parametric Mapping software (SPM8, http://www.fil.ion.ucl.ac.uk/spm) 



The preprocessing included slice-timing correction relative to middle ax- 
ial slice for the temporal difference in acquisition among different slices, 
head motion correction, spatial normalization into the Montreal Neurolog- 
ical Institute stereotaxic space, resampling to 3-mm isotropic voxels. 8(9) 
subjects were excluded from the dataset with TR=0.645s (TR=1.4s) be- 
cause either translation or rotation exceeded ±1.5 mm or ±1.5°, resulting in 
16(TR=0.645s) and 15(TR=1.4s) subjects each one scanned in two sessions 
which were used in the analysis). One subject whose data were too noisy was 
excluded from the TR=3 dataset, resulting in 10 subjects used in the analy- 
sis. In order to avoid introducing artificial local spatial correlations between 
voxels, no spatial smoo t hing was applied for further analysis, as pre viously 



suggested (jZhang et all 1201 It I Salvador et al.l . 120051 ; iLiao et al.l . l201ll ) . 



2.6. Anatomical parcellation and analysis 

The functional images were segmented into 90 regions of interest (ROI) 
using automated anatomical labeling (AAL) template as reported in previous 
studies. For each subject, the representative time series of each ROI was 
obtained by ay e r aging the fMRI time series across all voxels in the ROI 



([Salvador et al.l . 120051 ). Several procedures were used to remove possible 



spurious variances from the data through linear regression. These were i) 
six head motion parameters obtained in the realigning step, ii) signal from a 
region in cerebrospinal fluid, iii) signal from a region centered in the white 
matter, iv) global signal averaged over the whole brain. The BOLD time 



3 http://fcon_1000.projects. nitrc.org/indi/pro/eNKI_RS_TRT/FrontPage.html. 
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series were deconvolved into neural state signal using the above mentioned 
approach. 



2. 7. Effective connectivity network analysis 

The topological properties of the effective connectivity network were de- 
fined on the basis of a 90 x 90 binary directed graph G, consisting of nodes 
and directed edges: 



Where refers to the directed edge from ROI % to ROI j in the graph. T 
indicates the threshold. In a directed graph is not necessarily equal to e^. 
Considering that the graph we focused on is directed, all topological proper- 
ties were calculated on incoming and outgoing matrix, respectively. Graph 
theoretical analyses were carried out on the effective connectivity network 
using the Brain Connectivity Toolbox ( jRubinov and Spornsl . 120101 1 . 



2.8. Threshold selection 

As previous studies suggested that the brain networks of each subject 



normally differ in both the number and weighting of the edges ( iZhang et al. 



20 111 ; iLiao et al.l . l201ll ). we applied a matching strategy to characterize the 
properties of effective connectivity network. Both the global and local net- 
work efficiencies have a propensity for b eing higher with greater numbers 



of edges in the graph (IWen et all l201ll ). Modifying the sparsity values 



(number of edges) of the adjacency matrix also altered the graph's struc- 
ture. As a consequence it was suggested that the graphs to be compared 
must have (a) the same num ber of nodes and (b) the same number of edges 
(IBullmore and Bassettl . 1201 ll ) . The cost was defined as the ratio of the num- 
ber of existing edges divided by the maximum possible number of edges in 
a network. Since there is currently no formal consensus regarding selection 
of cost thresholds, here we selected a range of 0.05 to 0.14 with step = 0.01 
for subsequent network analyses. The lower bound was chosen as the one 
yielding a sparse graph with mean degree > 2Zn(90) (total number of edges 
> 405 where 405/90 2 = 0.05). The upper threshold corresponded to the 
smallest significant value of Granger causality (F-test with p = 0.05) across 
all subjects). 
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2.9. Network metrics 

For effective connectivity network at each cost threshol d, we calculated 
both overall topological properties and nodal characteristics ( IRubinov and Sporns . 



2010l ). The overall topological properties included i) small- worldness (a), re- 



lated to normalized clustering coefficient (7) and normalized characteristic 
path length (A); ii) network efficiency, divided in local efficiency (Ei oc ) and 
global efficiency (E g i b). The nodal characteristics included i) the nodal de- 
gree, that quantifies the extent to which a node is relevant to the graph, and 
ii) the nodal efficiency, that quantifi es the importance of the no des for the 
communication within the network ( iBassett and Bullmord . 120061 ) . Further- 
more we calculated the area under the curve (AUC) across all cost values for 
the above mentioned network properties. This quantity represents a summa- 
rized scalar for topological characterization of brain networks independent of 
single cost threshold selection. 



3. Results 



3.1. Reconstruction of HRF 

We tested the proposed deconvolution method on resting-state fMRI data; 
following the procedure summarized in the box, firstly we set a maximum 
time lag from a given threshold crossing, and obtain an optimal value for this 
lag, denoted with k. The histograms for k, reported in Fig J2] show a maximum 
around 4 ~ 6s , which is consistent with a p revious study acc ording to which 



the latency delay is 4 ~ 8s in gray matter ( ILee et al.l . Il995h . It is worth to 



mention that the lower TR could allow a more accurate estimation of the lag. 



To assess the effect of deconvolution, we compared the shape of voxel 
based HRF over the whole brain using different TRs. We focused on three pa- 
rameters: response height, time-to-peak, and full-width at half-max (FWHM) 
as potential measures of response magnitude, latency, and duration. Us- 
ing principal component analysis we determined the average intersubject 
variability of HRF maps. We found that the first component of HRF ac- 
counted for 81.7 ± 2.9%(response height), 98.1 ± 1.2%(time to peak) and 
95.6 ± 3.5%(FWHM) of the variance. Furthermore, the spatial distribution 
is very similar to the mean group map. The mean group results are plotted in 
Figj3j The response height, time to peak and FWHM of HRFs differ across 
brain regions, as a consequence of multiple factors including neural activ- 
ity differences, global magnetic susceptibilities, vascular differences, baseline 
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Figure 2: Histogram of the values of time deviation k. (*:without regression of global 
signal) 



cerebral blood flow, slice timing differences etc. (IHandwerker et al.l . 120041 ) . 
These patterns are remarkably similar across subjects and TRs, reflecting 
the robustness of the proposed approach. 



3.2. Variance stability of causality matrix 

As another indicator of the stability of the proposed joint of deconvolu- 
tion and PCGC approach we tested the variance of causality matrix across 
all subjects. We calculated the variance of the Granger causality matrix ob- 
tained both on the BOLD and deconvolved BOLD level signal. Firstly, we 
converted the matrix to Z-scores, then we calculated the variance of each 
matrix element, finally summing up the all these values into an overall vari- 
ance index. The variance of Granger causality matrix obtained from the 
deconvolved signal is much lower than the one of the BOLD level matrix for 
all TR values (Figj4j). Also, PCGC method kept the variance lower than 
full conditioned GC method. This result was confirmed testing a network 
at another scale using 1024 nodes (FigJU the native AAL segmentation was 
parcellated into 1024 micr o regions of interes t of approximately identical size 



across both hemispheres ( jZhang et al.l . 1201 ll ); in this case we could not test 



full conditioned GC due to small number of samples). 

3.3. Global signal regression 

As shown in previous studies, several sources of spurious variance should 
be removed by regression: motion artifacts, white matter and ventricular 
time courses. Still, the effects of regression against the global signal, calcu- 
lated by averaging across all voxels within a whole brain mask, are debated. 
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Figure 3: Spatial distribution of HRF shape: response height, time to peak and full- 
width at half-max (all values have been normalized, keeping range from to l).(*:without 
regression of global signal) 
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Figure 4: Total variance of causality matrix across all subjects. Full conditional Granger 
causality (CGC) and PCGC combined with BOLD and deconvolved BOLD level signal 
were used for construction of causality matrix. 



In order to evaluate this effect on our data we calculated spatial correla- 
tion between the the group mean image of HRF (response height, time to 
peak, FWHM) with and without regression of global signal in the preprocess- 
ing step, obtaining high Pearson correlation between them: r=0.97(response 
height), 0.90(time to peak), 0.88(FWHM). We can thus conclude that re- 
gression against global signal still preserved the spatial distribution. 

3.4- Effective connectivity network recovery with partial conditioning 

When trying to reconstruct effective connectivity networks, we are faced 
with the problem of coping with a large number of variables, when the appli- 
cation of multivariate Granger causality may be questionable or even unfeasi- 
ble, whilst bivariate Granger causality would detect also indirect interactions. 
Conditioning on a large number of variables requires an high number of sam- 
ples in order to get reliable results. Reducing the number of variables that one 
has to condition over would thus provide better results for small data-sets. In 
the general formulation of Granger causality, one has no way to choose this 
reduced set of variables; on the other hand, in the framework of information 
theory, it is possible to individuate the most informative variables one by one. 

The optimal model order m (order of the autoregressive model in Granger 
causality, embedding dimension in transfer entropy) for deconvolved BOLD 
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and BOLD signal was determined by leave-one-out cross-validation, and was 
found to be 3 for TR=0.645s, 2 for TR=1.4s and 1 for TR=3s. Under the 
Gaussian assumption, we constructed effective connectivity network using 
PCGC method. We firstly have to determine the number of variables upon 
which conditioning. To do this we look at how much uncertainty we elimi- 
nate adding an extra variable, letting the number of conditioning variables 
included rid vary from 1 to 20. This uncertainty can be expressed in terms 
of the information that we gain adding an extra variable. In Figj5l we plot 
the information gain as a function of rid] as expected, both this quantity and 
its increment decrease monotonically with rid- 



We can observe that the knee of the curves occurs when six variables 
are considered. This happens also when we consider different brain prior 
templates with 17 or 160 nodes (results not reported here). This could be 
connected to the fact that the average number of modul es which explain 



the equal-time correlations of resting brain is close to six ( iMarinazzo et al. 



20101 ; ISalvador et all 120051 ). therefore picking one variable from each module 



is sufficient to have most of the information, about a given channel, that 
can be obtained from the remaining channels, and this independently on the 
number of nodes. The effect of deconvolution, for the information gain, is to 
qualitatively raise the curve for TR=0.645s, and to lower them for TR=1.4s. 
This trend (not statistically significant) might be the result of two competing 
effects, the fact the deconvolution may remove spurious correlations and/or 
restore genuine correlations obscured by noise. 



Synthesizing the knee of the curves, sensitivity and specificity, we con- 
sider rid = 10 as the most appropriate number of most informative variables 
to include in the conditioning procedure. 



3.5. Global characteristics of ECN 

The global topological properties of brain ECN at deconvolved BOLD and 
BOLD level rely on the choice of thresholds. We used multiple cost thresholds 
and the AUC to evaluate the stability of the topological organization (Table 
1). An higher number of differences between the two networks was found 
with a (relatively) longer TR (TR=1.4s). Specifically, the AUC of small- 
worldness (a), normalized clustering coefficient (7), clustering coefficient (C p ) 
and local efficiency(£'; oc ) displayed the most significant differences, similar to 
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Figure 5: The mutual information gain (Ay), when the (rid + l)-th variable is included, is 
plotted versus rid- The information gain is averaged over all the variables. 
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Global network parameter 


AUC difference 


TR=0.645s 


TR=1.4s 


Incoming 


Outgoing 


Incoming 


Outgoing 


Sigma(cr) 


Y 


N 


Y 


Y 


Lambda(A) 


N 


N 


N 


N 


Gamma(7) 


Y 


N 


Y 


Y 


Characteristic path length(L p ) 


N 


Y 


Y 


N 


Clustering coefficient (C p ) 


Y 


N 


Y 


Y 


Global efficiency (E g i ob ) 


N 


Y 


Y 


N 


Local efficiency(i^ oc ) 


Y 


N 


Y 


Y 



Table 1: Comparison of AUC between deconvolved BOLD and BOLD. n c i = 10, Y: p < 
0.05, FDR corrected; N: otherwise. 



what emerged with TR=0.645s. For the data set with shorter TR we found 
significant differences in the characteristic path length and global efficiency 
of the outgoing network, whereas the most relevant differences were found 
for the incoming network with the longer TR. 

3. 6. Nodal characteristics of ECN 

Comparing the two networks on nodal degree, nodal global efficiency 
and nodal local efficiency revealed modifications in deconvolved BOLD and 
BOLD level (Figj6]). The patterns of nodal degree modifications resembled to 
those of nodal global efficiency in incoming network in both TR=0.645s and 
TR=1.4s fMRI data sets. In addition, more brain regions showed modified 
nodal degree and (global/local) efficiency in TR=0.645s rather than TR=1.4s 
data. 

4. Discussion 

In this study we proposed a novel methodology to achieve deconvolution 
in resting state data using spontaneous pseudo events, and to apply partially 
conditioned Granger Causality to the analysis of fMRI data. In our opinion 
this joint approach is the most convenient to infer effective connectivity with 
Granger Causality from resting state fMRI data. 

In the absence of a well defined ground truth, and in the light of the 
still active and unresolved debate on the usefulness of HRF deconvolution 
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Figure 6: Z-scores for Area under Curve from regional nodal parameters (deconvolved 
BOLD vs BOLD), p < 0.05, FDR corrected, (n<j = 10). Blue indicate negative values, red 
positive values. The point size is proportional to the absolute Z value. 



Granger causality based connectivity, we limit ourselves to validate the stabil- 
ity of the proposed method and indicate a possible path for the continuation 
of this debate, quantifying and comparing the overall topological properties 
of large-scale ECNs on deconvolved BOLD-level versus BOLD-level signals, 
investigating also the effect of different time resolutions (TR=0.645s and 
TR=1.4s). 



Previous discussions on evaluating effective connectivity from fMRI data 
reached the conclusion that it is bett er to use state-spa c e model for infer- 



ring c ausality on hidden neural states (IValdes-Sosa et all 1201 It iRyali et al. 



201 ll ; iBakhtiari and Hossein-Zadehl . l2012|). A pioneering EEG-fMRI study 
provided the first experimental substantiation of the theoretical possibility to 
i mprove interregion al coupling estimation from hi d den n eural states of fMRI 
( iDavid et all 120081 ) . Though promising (IFristonl . 120091 ) . these implications 
are still limited by the fact that multimodal recording is invasive and not 
applicable to healthy controls. As a consequence, data-driven methods for 
substantiating the confounding variability of haemodynamics have been de- 
veloped. The two availab le types of state space models in estimation of HRF 
(IValdes-Sosa et all 120111): t he generic (linear canonical/spline HRF) (jGloverl . 



19991: Marrelec et al 



20031) and biophysically informed models (DCM non- 
linear HRF) (IFriston et all |2000| ). Generic models are widely applicable but 
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lack specific biophysical constraints ( iGloverl . Il999t iMarrelec et al.l . 120031 ) . 
whil e biophysically infor med models are constrained by the hypothesis it- 
self (IFriston et al.l . l2000h . A recently proposed, biophysically informed bind 
deconvolution approach based on the state-of-the-art Cubature Kalman fi l- 



tering could be a useful tool for resting-state fMRI (IHavlicek et al.l . 1201 ll ). 
In the present study, however, we use a simpler approach which employs 
the generic linear canonical HRF for deconvolution. It is worth to point 
out that the significant differences between BOLD- and deconvolved BOLD- 
level effective connectivity found in complex network measures cannot abso- 
lutely exclude the misestimation of HRF. Furthermore HRF latency effect 
does not always critically affect the evaluation of mutual influence, so ECNs 
on BOLD and deconvolved BOLD level could have important consistencies 
flSupekar and Menonl . lioij ). 



Findings from brain connectivity studies have now demonstrated that the 
human brain network exhibits robust small-world topological properties, not 
only in the anatomical connectivity (reconstructed by diffusio n tractography) 



dHagmann et al.l . 120081 ) and functional connectivity network ([Salvador et al 



20051 ). but also in effective connectivity network (ILiao et al.l . l201lh . The 



current results also suggested that the ECNs obtained from BOLD and de- 
convolved data, with shorter and longer TR, have prominent small-world 
attributes, which would thus be confirmed as a general signature of robust 
organization of complex brain networks. Small-worldness indicates indeed 
an optimal balance b etween segregated and integ rated organization to pro- 
cess the information ( iBassett and Bullmord . 120061 ). For relatively longer TR 
we found significant differences between BOLD and deconvolved ECNs. Al- 
though an explanation based on precise neurobiological mechanisms is still 
not evident, we can suggest that the BOLD effect results from a more com- 
plex se quence of effects l inking neuronal activity, vascular changes and MRI 
signal (iLogothetisl . 120081 ). Hemodynamic delay, and hen ce the corre c t ons et 



of the events is indeed hard to capture with a long TR (ILaufs et al.l . 120081 ). 



In complex networks organization, the normalized clustering coefficient 
and the clustering coefficient are two key measures. They quantify the ex- 
tent of l ocal cliquishness or of local efficiency of information transfer of a 
network ( iBullmore and Bassettl . 1201 ll ). reflecting the local properties of net- 
work topologies. For longer TR, we observed significant differences between 
the two level ECNs. Thus the short-scale or local-scale network proper- 
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ties are indeed affected by deconvolution. Moreover, the normalized char- 
acteristic path length and the characteristic path length quantify global 
efficie ncy or the capability f or pa rallel information propagation of a net- 
work (jBullmore and Spornsl . 120091 ). These two measurements along with 
global efficiency are mainly associated with long-range co nnections ensur- 



ing effective interactions or rapid transfers of information (]He et al.l . 120091 ) 



It is widely accepted that long-range axonal connectivity being an impor- 
tant indj£atoj^_of^he^^ organization of the human cor- 
tex (IKnosche and Tittgemeyerl . l201ll ). This study reported no differences in 
long-range network organization. 



It is known that resting-state functional connectivity studies using ei- 
ther seed functional connectivity or independent component analysis bene- 
fit from higher sam pling rates to adeq uately sample undesirable respiration 



and cardiac effects (IBirn et all 120081 ) . while for event-related fMRI, faster 



sampling could allow for a better characterization of the hemodynamic re- 
sponse. The same applies to GCA. The previous simulations showed that 
accuracy of Granger causality depends on volume TR, faster sampling inter- 
yal increased the detection capacity of GCA of fM RI data to neural causality 
( iDeshpande et all . |2010| ; iRoebroeck et al.l . 120051 ). In this paper, we focus on 
resting-state fMRI data with TR=0.645s and 1.4s to maximally escape infor- 
mation loss due to low sampling. Considering the limitation of acquisition 
sequence, the conventi onal fast TR data acquisition br i ngs t o the loss of the 



fine spatial resolution (jHuettel et al.l . l2004t iKim et al.l . |l994J) . 



Other methodological considerations are worth to be mentioned. The 
first one concerns data preprocessing. As a general idea spatial smoothing 
can reduce the noise and increase signal-t o-noise ratio, t herefo re improving 
the accuracy of detecting of neural event (jHuettel et al.l . 120041 ) . Here we do 



not include this step. As we used AAL template, spatial smoothing would 
blur the boundary among these regions, which may affect the GC inference. 
Temporal filtering is frequently a necessary step for functional connectivity 
analysis of resting-state fMRI data . In line with previou s studies that con - 



sidered a low model order in GCA ( [Hamilton et al.l . |2010| ; iLiao et al.l . 1201 ll ) 



we did not performed low-pass filtering. 



Secondly, graph theoretic approach is one of the most powerful and flex- 
ible approaches to investigate functional and structural brain connectome; 
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still some controversies remain, concerning; the definitio n of network nodes 
and edges (IBullmore and Bassettl . l201ll ; IWig et al.l.l201ll ) . Diff erent node def- 
i nitions by prior ana t omic brain template s ( jWang et al.l . 120091 ) or node scales 



(jFornito et al.l . |2010| ; IZalesky et al.l . 120101 ) could produce different results. In 



future works, more brain templates and more node scales comparison for ef- 
fective connectivity network should be explored. 



Acknowledgments 

H. Chen was supported by the Natural Science Foundation of China (No. 
61125304 and No. 61035006). G.R. Wu gratefully acknowledges the financial 
support from China Scholarship Council(2011607033). 

References 
References 

Angelini, L., De Tommaso, M., Marinazzo, D., Nitti, L., Pellicoro, M., Stra- 
maglia, S., 2010. Redundant variables and Granger causality. Physical 
Review E 81, 037201. 

Bakhtiari, S., Hossein-Zadeh, C, 2012. Subspace-based identification algo- 
rithm for characterizing causal networks in resting brain. Neuro Image 60, 
1236 - 1249. 

Barnett, L., Barrett, A., Seth, A., 2009. Granger causality and transfer 
entropy are equivalent for gaussian variables. Physical Review Letters 
103, 238701. 

Barnett, L., Seth, A.K., 2011. Behaviour of Granger causality under filtering: 
theoretical invariance and practical application. Journal of Neuroscience 
Methods 201, 404 - 419. 

Bassett, D., Bullmore, E., 2006. Small-world brain networks. The neurosci- 
entist 12, 512-523. 

Birn, R., Murphy, K., Bandettini, P., 2008. The effect of respiration varia- 
tions on independent component analysis results of resting state functional 
connectivity. Human brain mapping 29, 740-750. 



21 



Bressler, S., Seth, A., 2011. Wiener-Granger causality: A well established 
methodology. Neuroimage 58, 323-329. 

Biichel, C, Friston, K., 1997. Modulation of connectivity in visual path- 
ways by attention: cortical interactions evaluated with structural equation 
modelling and fMRI. Cerebral Cortex 7, 768-778. 

Bullmore, E., Bassett, D., 2011. Brain graphs: graphical models of the human 
brain connectome. Annual Review of Clinical Psychology 7, 113-140. 

Bullmore, E., Sporns, O., 2009. Complex brain networks: graph theoretical 
analysis of structural and functional systems. Nature Reviews Neuroscience 
10, 186-198. 

Buxton, R., Wong, E., Frank, L., 1998. Dynamics of blood flow and oxy- 
genation changes during brain activation: the balloon model. Magnetic 
Resonance in Medicine 39, 855-864. 

David, O., Guillemain, I., Saillet, S., Reyt, S., Deransart, C, Segebarth, C, 
Depaulis, A., 2008. Identifying neural drivers with functional MRI: an 
electrophysiological validation. PLoS Biology 6, e315. 

Deco, C, Jirsa, V., 2012. Ongoing cortical activity at rest: criticality, multi- 
stability, and ghost attractors. The Journal of Neuroscience 32, 3366-3375. 

Deshpande, C, Sathian, K., Hu, X., 2010. Effect of hemodynamic variability 
on Granger causality analysis of fMRI. Neuroimage 52, 884-896. 

Feinberg, D., Moeller, S., Smith, S., Auerbach, E., Ramanna, S., Glasser, M., 
Miller, K., Ugurbil, K., Yacoub, E., 2010. Multiplexed echo planar imaging 
for sub-second whole brain fMRI and fast diffusion imaging. PLoS One 5, 
el5710. 

Fornito, A., Zalesky, A., Bullmore, E., 2010. Network scaling effects in graph 
analytic studies of human resting-state fMRI data. Frontiers in Systems 
Neuroscience 4. 

Friston, K., 2009. Causal modelling and brain connectivity in functional 
magnetic resonance imaging. PLoS Biology 7, el 000033. 

Friston, K., Harrison, L., Penny, W., 2003. Dynamic causal modelling. Neu- 
roimage 19, 1273-1302. 



22 



Friston, K., Mechelli, A., Turner, R., Price, C, 2000. Nonlinear responses 
in fMRI: the balloon model, volterra kernels, and other hemodynamics. 
Neurolmage 12, 466-477. 

Friston, K., Trujillo-Barreto, N., Daunizeau, J., 2008. Dem: a variational 
treatment of dynamic systems. Neuroimage 41, 849-885. 

Glover, G., 1999. Deconvolution of impulse response in event-related bold 
fMRI. Neurolmage 9, 416-429. 

Hagmann, P., Cammoun, L., Gigandet, X., Meuli, R., Honey, C, Wedeen, V., 
Sporns, O., 2008. Mapping the structural core of human cerebral cortex. 
PLoS Biology 6, el59. 

Hamilton, J., Chen, G., Thomason, M., Schwartz, M., Gotlib, I., 2010. Inves- 
tigating neural primacy in major depressive disorder: multivariate Granger 
causality analysis of resting-state fMRI time-series data. Molecular Psy- 
chiatry 16, 763-772. 

Handwerker, D., Gonzalez-Castillo, J., D'Esposito, M., Bandettini, P., 2012. 
The continuing challenge of understanding and modeling hemodynamic 
variation in fMRI. Neurolmage 62, 1017-1023. 

Handwerker, D., Ollinger, J., D'Esposito, M., 2004. Variation of bold hemo- 
dynamic responses across subjects and brain regions and their effects on 
statistical analyses. Neuroimage 21, 1639-1651. 

Havlicek, M., Friston, K., Jan, J., Brazdil, M., Calhoun, V., 2011. Dynamic 
modeling of neuronal responses in fMRI using cubature Kalman filtering. 
Neurolmage 56, 2109 - 2128. 

Havlicek, M., Jan, J., Brazdil, M., Calhoun, V., 2010. Dynamic Granger 
causality based on Kalman filter for evaluation of functional network con- 
nectivity in fMRI data. Neurolmage 53, 65-77. 

He, Y., Dagher, A., Chen, Z., Charil, A., Zijdenbos, A., Worsley, K., Evans, 
A., 2009. Impaired small-world efficiency in structural cortical networks 
in multiple sclerosis associated with white matter lesion load. Brain 132, 
3366-3379. 



23 



Huettel, S., Song, A., McCarthy, G., 2004. Functional magnetic resonance 
imaging, volume 1. Sinauer Associates Sunderland, MA. 

Kim, S., Hendrich, K., Hu, X., Merkle, H., Ugurbil, K., 1994. Potential 
pitfalls of functional MRI using conventional gradient-recalled echo tech- 
niques. NMR in Biomedicine 7, 69-74. 

Knosche, T., Tittgemeyer, M., 2011. The role of long-range connectivity 
for the characterization of the functional-anatomical organization of the 
cortex. Frontiers in systems neuroscience 5. 

Laufs, H., Daunizeau, J., Carmichael, D., Kleinschmidt, A., 2008. Recent ad- 
vances in recording electrophysiological data simultaneously with magnetic 
resonance imaging. Neuroimage 40, 515-528. 

Lee, A., Glover, G., Meyer, C, 1995. Discrimination of large venous vessels in 
time-course spiral blood-oxygen-level-dependent magnetic-resonance func- 
tional neuroimaging. Magnetic Resonance in Medicine 33, 745-754. 

Liao, W., Ding, J., Marinazzo, D., Xu, Q., Wang, Z., Yuan, C, Zhang, Z., 
Lu, G., Chen, H., 2011. Small-world directed networks in the human brain: 
multivariate Granger causality analysis of resting-state fMRI. Neuroimage 
54, 2683-2694. 

Logothetis, N., 2008. What we can do and what we cannot do with fMRI. 
Nature 453, 869-878. 

Marinazzo, D., Liao, W., Pellicoro, M., Stramaglia, S., 2010. Grouping time 
series by pairwise measures of redundancy. Physics Letters A 374, 4040- 
4044. 

Marinazzo, D., Pellicoro, M., Stramaglia, S., 2012. Causal information ap- 
proach to partial conditioning in multivariate data sets. Computational 
and Mathematical Methods in Medicine 2012, 8. 

Marrelec, G., Benali, H., Ciuciu, P., Pelegrini-Issac, M., Poline, J., 2003. Ro- 
bust bayesian estimation of the hemodynamic response function in event- 
related bold fMRI using basic physiological information. Human Brain 
Mapping 19, 1-17. 



24 



Petridou, N., Gaudes, C, Dryden, I., Francis, S., Gowland, P., 2012. Periods 
of rest in fMRI contain individual spontaneous events which are related to 
slowly fluctuating spontaneous activity. Human Brain Mapping, in press 
(DOI: 10.1002/hbm.21513). 

Riera, J., Watanabe, J., Kazuki, I., Naoki, M., Aubert, E., Ozaki, T., 
Kawashima, R., 2004. A state-space model of the hemodynamic approach: 
nonlinear filtering of bold signals. Neurolmage 21, 547-567. 

Roebroeck, A., Formisano, E., Goebel, R., 2005. Mapping directed influence 
over the brain using Granger causality and fMRI. Neuroimage 25, 230-242. 

Roebroeck, A., Formisano, E., Goebel, R., 2011. The identification of inter- 
acting networks in the brain using fMRI: model selection, causality and 
deconvolution. Neuroimage 58, 296 - 302. 

Rubinov, M., Sporns, O., 2010. Complex network measures of brain connec- 
tivity: uses and interpretations. Neuroimage 52, 1059-1069. 

Ryali, S., Supekar, K., Chen, T., Menon, V., 2011. Multivariate dynamical 
systems models for estimating causal interactions in fMRI. Neuroimage 
54, 807-823. 

Salvador, R., Suckling, J., Coleman, M., Pickard, J., Menon, D., Bullmore, 
E., 2005. Neurophysiological architecture of functional magnetic resonance 
images of human brain. Cerebral Cortex 15, 1332-1342. 

Smith, S., Miller, K., Salimi-Khorshidi, G., Webster, M., Beckmann, C, 
Nichols, T., Ramsey, J., Woolrich, M., 2011. Network modelling methods 
for fMRI. Neuroimage 54, 875-891. 

Stephan, K.E., Roebroeck, A., 2012. A short history of causal modeling of 
fMRI data. Neuroimage 62, 856-863. 

Supekar, K., Menon, V., 2012. Developmental maturation of dynamic causal 
control signals in higher-order cognition: A neuro cognitive network model. 
PLoS Computational Biology 8, el002374. 

Tagliazucchi, E., Balenzuela, P., Fraiman, D., Chialvo, D., 2012. Criticality in 
large-scale brain fMRI dynamics unveiled by a novel point process analysis. 
Frontiers in Physiology 3. 



25 



Valdes-Sosa, P., Roebroeck, A., Daunizeau, J., Friston, K., 2011. Effective 
connectivity: Inffuence, causality and biophysical modeling. Neuroimage 
58, 339 - 361. 

Wang, J., Wang, L., Zang, Y., Yang, H., Tang, H., Gong, Q., Chen, Z., 
Zhu, C, He, Y., 2009. Parcellation-dependent small-world brain functional 
networks: A resting-state fMRI study. Human Brain Mapping 30, 1511- 
1523. 

Wen, W., Zhu, W., He, Y., Kochan, N., Reppermund, S., Slavin, M., Brodaty, 
H., Crawford, J., Xia, A., Sachdev, P., 2011. Discrete neuroanatomical 
networks are associated with specific cognitive abilities in old age. The 
Journal of Neuroscience 31, 1204-1212. 

Wig, C, Schlaggar, B., Petersen, S., 2011. Concepts and principles in the 
analysis of brain networks. Annals of the New York Academy of Sciences 
1224, 126-146. 

Zalesky, A., Fornito, A., Harding, I., Cocchi, L., Yiicel, M., Pantelis, C, 
Bullmore, E., 2010. Whole-brain anatomical networks: Does the choice of 
nodes matter? Neuroimage 50, 970-983. 

Zhang, Z., Liao, W., Chen, H., Mantini, D., Ding, J., Xu, Q., Wang, Z., 
Yuan, C, Chen, C, Jiao, Q., et al., 2011. Altered functional-structural 
coupling of large-scale brain networks in idiopathic generalized epilepsy. 
Brain 134, 2912-2928. 



26 



